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Abstract. Dynamical evolution of stellar mass distribution in star clus¬ 
ters is analysed by considering simultaneously the effects of dynamical friction, 
stochastic heating and the gravitational potential due to mass distribution in 
the clusters. A simple expression is suggested for the dynamical friction which 
adequately describes it in both high and low velocity ranges and the effect of 
mass distribution in the cluster on the dynamics of a test mass can be described 
by an anharmonic potential. With the help of energy considerations we describe 
the dispersions in position and in velocity for stars of different masses show¬ 
ing how the slowing down and mass segregation of stars evolve in the cluster. 
Results are presented with the parameters for the star clusters of our Galaxy. 
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1. Introduction 

Properties of star clusters form some of the most fascinating areas of astrophysics, partic¬ 
ularly due to their close contact with issues related to star formation, mass function and 
stellar evolution. Such systems have stars over a wide range of masses (m) which lie in the 
range 0.1 < m/MQ < 10. The distribution of these masses in real space as also in velocity 
space are signatures of various effects that arise in the context of dynamical evolution of 
the star in the gravitational potential of all the masses in the cluster. In what follows, 
we give a simple quantitative estimate of the evolution of velocity of any star of mass M, 
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(for the moment only those whose life times are longer than the age of the cluster). The 
method is completely analytical and thus amenable to simple computational schemes. 
This is made possible firstly because a simple and yet accurate approximations for the 
Chandrasekhar dynamical friction formula could be found and the same type of formulae 
also found to fit the gravitational field in the cluster. This choice enables us to integrate 
the equation of motion analytically and its usefulness with reference to stellar dynamics 
in a spherically symmetric potential is demonstrated in the subsequent sections with spe¬ 
cial reference to open clusters and globular clusters. The relaxation process is considered 
to be due to dynamical friction, the violent relaxation process being completed within 
few thousands of years. The mass burning upto this time since the birth of the star is 
indeed negligible and u(0), M(0) are thus the quantities that obtain after the completion 
of violent relaxation. The gas and dust content in clusters being small, the drag forces 
due to dust and gas are neglected and only dynamical friction due to stars is considered. 
As is known, the dynamical friction is caused by the fluctuations in the trajectory of the 
star in question, due to encounters (not physical collisions) with other stars present in 
the cluster. The other important parameter that contributes to the dynamics of a test 
particle is the gravitational held arising out of the mass distribution (assumed Gaussian) 
in the cluster. We show that for a density profile p(r) decaying with r, the gravitational 
potential ())(r) too can be adequately described by a function similar to that used to ap¬ 
proximate the Chandrasekhar dynamical friction expression. A combination of these two 
approximations helps to integrate the equation of motion through simple mathematical 
operations. The integrated equation thus considers the dynamical evolution of a test mass 
M due to combined effects of dynamical friction and its oscillation in an anharmonic 4>{r) 
due to mass distribution in the cluster. It is indeed necessary to explain the motivation 
behind the present analytical work. This is mainly motivated by the recognition that 
’’there is plenty of room for the development of new analytic techniques”, ’’bringing into 
play more physically motivated reasoning in the choice of fitting functions” (Hut, 1996). 
The results that follow do give tractable analytical methods to estimate various physical 
processes in the evolution of a test mass in a cluster and may prove to be useful for the 
future theoretical understanding of the existing observational data. 

The main principles involved in the dynamics of the system are discussed in the next 
two sections. In section 2, we determine the effect of dynamical friction by considering 
stellar encounters. Essentially, the method is similar to what is given by Binney and 
Tremaine (1987) in which a King type velocity distribution has been incorporated. The 
main scheme for integrating the equation of motion is contained in sections 3 and 4 
while the numerical scheme of solution is given in section 5. An interesting feature of 
our scheme is that it can easily accommodate modification due to mass loss arising out 
of various mechanisms that appear in the course of stellar evolution and is described in 
section 6. The effect of the above phenomena in describing mass segregation in the cluster 
is described in section 7 which is followed by results and discussions in section 8 of the 
paper. 
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2. Dynamical friction with King type velocity distribution 

The role of gravitational interaction between the masses is considered in terms of stellar 
encounters. This, as is understood from the celebrated work of Chandrasekhar leads to 
(1) slowing down of the velocity of the particle in the direction of its initial velocity uy (0) 
and (2) enhancements of its velocity vj_, in the component perpendicular to r||. The role 
of interparticle interaction is incorporated through the encounters between the test mass 
M, travelling with velocity vm and the scatterers m, travelling with velocity Vm- As is 
well known, the inverse square nature of the interaction leads to a change Avm, which 
is given by. 


^ 2m\Vo\cx r 

(M + to) [ G2(M + to)2 


( 1 ) 


where V^o = — vm and the components of unit vector a in the direction of initial Vq 

is cos 9q and in direction perpendicular to Vq is sin 9q, where cos 9q and sin 9q are given 
as G{m + M)/y/lG‘^{'m + M)2 + &2y4^ /y/(G'^{m + M)2 + b’^V^) respectively. 

These yield. 


|Avm||P 


(m + M)2 




|Aum_lP 


Am'^V^b'^ r &2^^4 

G2(m + M)4 [ G2(m + M)2 


which gives, 


( 2 ) 


( 2 . 1 ) 


\Avm\^ 


4mVo^ 
(m + M)2 



G2 (m + M) 



( 2 . 2 ) 


In order to determine the changes that occur in time t, we consider the encounter process 
to be uncorrelated and on an average 2TmoVobdb encounters occur per unit time for 
particles passing within an interval b and b + db. The physical situation puts the limits 
of b as bmin A b < b^ax and the velocity distribution of m follows a truncated Gaussian, 


f{Vm) 


f Ano 

1 Ir', 1 2 /o 2 ’ 

2a^ _ -v^ 2a^ 

1 0 



for V < V, 
for V > V, 


(3) 
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where Vg is the escape velocity of the scatterers from the cluster and A is a normalisation 
constant, being given by, 


iM 




(3.1) 


and we have defined x = VmlcrmV^ and Xg = VgfamV^- 

The velocity distribution of the scatterers as given in equations (3, 3.1) indeed satisfies 
the consistency condition that the cluster contains only those stars with velocity which 
cannot escape from the potential well created by the mass distribution in the cluster. The 
estimate for Vg can be easily made by integrating the Poisson equation and putting Ug/2 = 
(/)(oo), being thus consistent with the nature of the mass distribution. For an accurate 
determination of the dynamical friction one must keep in mind the mass spectrum and 
the mass dependence of the velocity dispersion am- In the scheme presented below, we 
assume that the mass spectrum can be represented by an average mass m and an average 
velocity dispersion am- 

It then follows from the equations (1), (2.2) and (3) that on account of encounters there 
appears a diminution of vy being given by, 


dVM 

dt 



AVm 2'Kbdb Vo /(wm) ■iTTVmdVm 


27rG^TO(m + M) 

2V2al 


[ln(Xmaa:/Xmm)] ^no(TO) VMFgh{x) 


( 4 ) 


d\AvM? 

dt 



\Avm\‘^ 2TTbdbVo f{Vm) 


\piam 


j X Fgfi{x^ 


(4.1) 


with X — ^ + + FlY and for all practical purposes we shall approximate 

ln(xmaa;/Xmm) « 2 In A with A = bmaxibmin- The Fgh{x) is given as 


(erf(a;) —-^e ~ “»)/a;^ ior x < Xg 

(erf(a;) — jx^ for x > Xg 


Fch{x) 


(4.2) 
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It can be seen that the above equations converge to the Chandrasekhar formulae when 
Xe ^ 1- In the low velocity range the dynamical friction oc vm and goes as vm/\vm\^ 
for the high velocity range. Further, it can be seen from simple algebra that v'j^ evolves 
as, 


vlfit + 6t) = vlt{t) + 2vM{t){dvM/dt)t5t + {d{vlj{t)) / dt)t5t (5) 


The second and third terms in equation (5) are consistent with the cooling term and 
heating term respectively in Bekenstein and Maoz (1992) ( see equations 3.2.5 and 4.16, 
therein) and are also consistent with equations (161-163) in Nelson and Tremaine (1999). 
In a future work we plan to re-examine their works in the light of equations (9-11) in the 
present paper. From equations (4.1) and (4.2) it can be seen that the second term in (5) 
is always negative while the last term in (5) is positive. The former gives a slowing down 
due to dynamical friction being oc m(rn + M) while the latter gives stochastic heating and 
varies as m^. Both originate in the inverse square interaction and act jointly, taking the 
system to equilibrium, being thus a manifestation of the fluctuation-dissipation theorem. 
For accurate evolution of the velocity we have to follow the formalism given in the next 
section. Full details of these calculations will be puplished shortly. We, however, present 
a summary in the next section in view of its importance on the question of evolution to 
equipartition. 


2.1 On approach to equipartition 

The question concerns, how the kinetic energy evolves in presence of encounters. 

It can be shown that on a single encounter the change in the kinetic energy of the mass 
M is given by (Binney and Tremaine 1987), 


A{Mv'l^/2) = (M/2) [2vm-Avm + \Avm\'^ + “^Avm-Vc] 


( 6 ) 


where vm is the change of velocity on encounter and Vc = {mVm + MvM)/{m + M) is 
the velocity of the center of mass of the system. Knowing the components Au|| and Ar;_L 
as given in equations (1-3) we Hnd the r.h.s. of equation (6) to be, 


-(M/2) 


(^'M|| -Wm||)(MuM|| +mv^ll) 


2G'^{m + M) 
G 2 (TO-kM) 2 -H 62 V^(/ 


(^M|| - 


2bG 

G^{m + M)^ + 


Vm\\)^{MvM± +mVjn±) 


(7) 
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If we have a situation in which particles m randomly encounter the test particle M at 
values of b which are +ve and —ve equally frequently, we have to use the averages (6) = 0 
and {vM\\Vm\\) = 0 (e.g. b +ve means that the projectile approaches from above and it’s 
track is bent downwards, while b —ve means that it approaches from below and the track 
is bent upwards, which essentially makes Aum_l to be +ve and —ve equally frequently, 
being thus zero on an average), so that 


^{Mvlt/2) 


\Mm/{m + M)^ 


„2 

-’Mil 


2G'^{m + Mf 

G‘^{m + My + b‘^V^ 


( 8 ) 


It shows that Mv]^/2 increases for and decreases for Mv'^n > so 

that in the former case the kinetic energy of M increases with scattering while it decreases 
in the latter case. Since there are 27r& db Vonf{vm) collisions per unit time, we find 


{d/dt){Mv1f/2) = —27r^G^mM(Mv^l|/i — ml 2 )n 


(9) 


where 


h = 


'o Vjf 


In 




[GHM + m)^ + bl,^Vo^\ 




( 10 ) 


72 = / AlA In 

Jo 


G^{M + mr+bl,,V,^ 

G2(M + m)2 + 6^.„y(4j 


/(vm)47ru^dv„ 


( 11 ) 


Dimensionally we note that I 2 /I 1 — Equation thus expresses the evolution of the 

kinetic energy Mv‘j^/2 of the test star and shows that its rate of change vanishes when the 
equipartition i.e. = m(u^l|) is reached. Full details of this evolutionary process, 

particularly, with respect to vm will be presented in a future publication. 

It is to be noted that the Chandrasekhar formula of which equations (4.2) and (13.2) 
are variants actually considers the density to be uniform. We have applied here the same 
formalism to an r dependent density profile in view of Kandrup’s result (Kandrup 1981) 
that in an r dependent density profile the Chandrasekhar formula can be adopted by 
using the local number density n(r) for the average n, appearing in the Chandrasekhar 
formula . Thus on averaging over all possible orbits (that lie at different r’ s) one may 
use the averaged number density, as we have done, to estimate the net effect of dynam¬ 
ical friction. Chandrasekhar’s simplification in considering the initial um( 0) and final 
vm{oo) to be straight lines have also been investigated and the corrections are found 
to be marginal (Kandrup 1983). Furthermore, Del Popolo and Gambera (1999) have 
shown that the numerical value of the dynamical friction can be changed from the Chan¬ 
drasekhar formula if a power law dependence n{r) ~ r~P (p > 0 see equation (12) in this 
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paper) be considered. We have , however, not incorporated this effect since such a power 
law dependence makes the number density n{r) blow up for r —> 0. In the computed 
results, given below, all the parameters have been taken from well known observations. 


3. Principles and Approximations 


We consider a test mass M, with an initial velocity v(0), undergoing motion in the 
cluster. This motion takes place in a gravitational field 4>{r,t) = 4>oir,t) + S<j){r,t) of 
which the first part describes slowly varying (in time) average gravitational field due to 
all stars in the cluster, while 6(j){r,t) is a rapidly varying time dependent potential that 
arises due to density fluctuations in the cluster owing to random velocities of stars i.e., 
d(j){r,t) gives rise to dynamical friction. The equation of motion of the star consists of 
the equations dL/dt = —L [Fdyn{M,v)/\v\] and dL^/dt = — [Fdyn{M,v)/\v\] for the 
angular momenta and f = —d(j)o/dr+{L'^/mr^) — VrFdyn{M, v)/|r^| for the radial motion. 
In what follows we consider L = 0 i.e. a star with no angular momentum, oscillating 
radially across the cluster. The case for non zero L will be dealt with separately in a 
future work. For the case at hand, the equation of motion (we consider only the radial 
velocity) thus reads, 

— = -VMr,t) - = F,ir)+6F,{r,t) - F^y^{M,v) (12) 


dv 

dt 


-^Mr) 


V 

t(M, v) 


( 12 . 1 ) 


where the first term in the r.h.s. of equation (12) gives the steady gravitational potential 
field, while the second term is a time dependent field arising due to changes in the 
density profile of the cluster in course of its dynamical evolution and the third term is the 
damping due to dynamical friction. By neglecting the SF^lr, t) we obtain equation (12.1). 
The field SFQ{r,t) though time dependent is a slowly varying quantity as compared to 
W6(j){r, t), which finally leads to dynamical friction. The net effect of the S(j){r,t) term 
has been employed here by introducing the last term i.e. dynamical friction in equations 
(12) and (12.1). We must note, however that the dynamical friction does have a role in 
the origin of the SFo{r,t) term also since the density profile evolves with time due to 
dynamical friction. Due to the slow variation of 6FQ{r,t), this term is ignored in this 
zeroth order approximation, to be, however, considered in a more complete theory. 

For integrating equation (12.1) two cases have to be considered separately. Equation 
(12.1) can be solved in its completeness by the method given by Van Kampen (1985). The 
method essentially involves the identification of a smallness parameter e and expansion of 
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the solution u as a power series in e. We find below that the system has a dimensionless 
parameter e' = relaxation time/time period of oscillation. When relaxation are slow we 
have e = 1/e' and the star behaves as a damped oscillator. For fast relaxation e" <C 1 and 
one can identify e = e' as the small parameter and the system behaves as an overdamped 
oscillator and hence does not demonstrate any to and fro motion. Considering the mass 
distribution in the cluster to be p{r) = poexp(—yr^) and a velocity distribution as given 
in equation (5) we find, 

-Fr = d(j)o/dr = yP{y) (13) 

where y = r^/y and 


1 

T{M,V,t) 


4TrG'^nom{m + M{0) — AM{t)) ln(A) 


Fch{x) 


(13.1) 


where Fch{x) has been defined in (4.4), AM{t) = \ f* {dM/dt) dt\ being the mass reduc¬ 
tion upto the time t after the birth of the star. It is clear that AM{t) depends upon the 
evolutionary phase of the star which is also determined by its initial mass M(0). In what 
follows we shall be concerned with the stars in the main sequence for which analytical 
expression for AM(0) will be used. It is seen that for a; ^ 0 Fch{x) ~ 4/3y^ while for 
X ^ oo Fch(x) ~ 1/a;^. We find that Fch{x) can be approximated as 


F{x) = l/[a + bx^]^^‘^ (13.2) 

with a = 3.25, b = 1.28 for globular cluster case and a = 2.56, b = 4.09 for open cluster 
case (see section 4) and , with a reliability of 99% as seen by Kolmogorov-Smirnov (KS) 
test, a simplification, which we shall use here for a very simple computational scheme. 
This approximation for F(x) also gives on integrating equation (13), and approximate 
expression for the gravitational potential inside the cluster, which reads 


<Po{r) = <PoHy) 

(13.3) 

where 


00 = 7r^/^Gpo/2y 

(13.4) 

II 

(13.5) 
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^(y) ~ y‘^F{y) + {3b/7)y'^[F{y)f + .... + .... 


(13.6) 


We were not able to find a closed form equation for (j){y) and the expansion given 
in equation (13.6) is fairly accurate upto y = 1.33, which will correspond to r ~ 40 pc 
for a typical globular cluster and r ~ 13 pc for an open cluster with halo. We note on 
numerically integrating (fiv) that it is monotonically increasing function and (fiy) —>-7/2 
for y oo. The escape velocity Ve of particles is then given by = (()o(^(oo) i.e. 
Ug = ^0 ■ In sll Ihs discussions that follow we are concerned with particles that are 
gravitationaly bound to the cluster and hence obey the restriction The choice 

of p(r) as a Gaussian is indeed arbitrary but it does satisfy that p(r) is peaked at r = 0 
and falls off rapidly beyond a length ro ~ 1/from the centre, which is an essential 
feature of the density profile of the cluster (Nilakshi et al. 2002). This results in a 
gravitational field Fr oc —r for r ^ tq and Fr oc —r~^ for r ^ rg. In what follows we 
ascribe such a behaviour through equations (13) and (13.2) The choices a and b that we 
present here are the ideal ones for the Gaussian p{r). However, for any other type of 
decay too similar approximations as given in equations (13) and (13.2) are possible in 
which the choices of a and b would be decided by the exact nature of the decay of p{r) 
e.g. for an exponential density distribution p{r) = p(0)e“”’’ we find that the Fo{r) can 
be approximated to Fq ~ {K,r)/[ai + foi(Kr)®]^/^ with 90% confidence as shown by the KS 
test. This observation thus demonstrates the usefulness of our approximation. For the 
subsequent calculation described here we continue to follow a Gaussian distribution for 
the density profile. 


4. Method of solution 

We note that the system introduces certain time scales. For x 0,F{x) 1, and for 

X oo,F{x) ~ (2/Y^)(2/3)a;^. Hence for low velocities we have a relaxational time 
scale t{M), 1/t{M,v) = (47rG^nom(m + M) ln(A)/2-\/2(T))^)(4/3-\/7r) and an oscillational 
time period Tq which is akin to the crossing time such that jT^ = (47r/3)Gpo- The 
solution of equation (12.1) must consider two distinct cases, (a) Tq <C t(M) and (&) 
To > r(M). 


4.1 Case (a) : Tq < t(M); e' = t{M)/Tq-, e = 1/e' 

We examine the problem in two cases. Firstly we consider a typical globular cluster with 
no = 25 pc“^, m = IMq, am = 6 km s“^, tq = 1/y^ ~ 30 pc with Ve = 22 km s“^. 
Next we consider the intermediate age open cluster NGC 2099 (M37) for which no « 3 
pc“^, m = IMq and am ~ 1 km s“^, ro = l/-y/7 ~ 10 pc, being the range upto which 
the halo is present. For this last case we find Ve = 2.6 km s“^ (Kalirai et al. 2001, 



10 


S. Chatterjee et al. 


Nilakshi and Sagar 2002). For stars with M(0) < IOMq one finds Treiax = 2r(M(0), 0) > 
9.2 X 10^ years while Tosc ~ 2.75 x 10^ years. Since Treiax/Tosc ~ 3, the star does exhibit 
oscillations, r{t) = a{t) cos{2TTt/Tosc) where r{t) is the instantaneous radial distance of 
the star from the center of the cluster. Indeed the nature of Fr shows the potential (j)Q{r) 
to be anharmonic, this anharmonicity can be incorporated in the calculation by making 
the oscillational time scales r(n)’s to be dependent on the amplitudes of the oscillation 
for the respective cycles n of oscillations of the particle. Indeed the complete dynamics 
requires the inclusion of higher harmonics of 2tt/T in) but as far as the energy loss is 
concerned it is the slowest oscillational mode that outweighs the contribution from other 
faster modes, due to rapid elimination of fast variables (Van Kampen 1985). It thus 
suffices to approximate the motion only in terms of the slowest 2tt/T in) mode as is done 
below. It is to be noted that with the passage of time the amplitude of oscillation a{t) 
will decay due to dynamical friction. Tosc is also a time dependent quantity arising due 
to anharmonicity of the potential 4>Q{r). If however, Treiax/Tosc < 1 the system becomes 
overdamped and such oscillations are not to be found as will be seen in case (b). 

The energy loss per cycle due to dynamical friction is then found by noting that power 
loss is the product of the velocity and the resistive force, so that on defining to be the 
energy of the particle in its n*^ cycle. 


En+l 


rT(n) 


— En — Ai?„ — 


= -2alX/f} 



FdyniM{t),v) V dt 

(m + MQ — AM t) cos^(27rt/T(n)) 
[1 + a;® cos®(27rt/T(n))]^/^ 


(14) 


where = (a/5)^/®, (3 = 47rG^noTOln(A)/2-\/2CT^/-ya, and x\ = u^/2ct^Vq, being 
the amplitude of the velocity in the nth cycle and T{n) is the corresponding time period. 
The above difference equation can be written as 




(^'n) — 


FoiZn) = 


= + + "o - ^m))T{n) 

,5 4 1965 4" 

32 [1 + z6]i/2 + 2097152 [1 + z6]i/2 

1_L TQarn \ 2 

10 ^ ” 


FNjZn) 

Fnizn) 


(15) 

(15.1) 

(15.2) 


where 4 = (48/15)^/®a;4 t = E = ( 7 / 2 ) = average radius of the cluster. 



Dynamics of star clusters 


11 


(27r/To)^ = (47r/3)Gpo) Tq being the oscillation time period of the star in the harmonic 
limit and t is the total time elapsed after the formation of the star. In equation (15) 
FN{zn) denotes diminution of dynamical friction as the star acquires larger velocity as it 
passes through the center of the cluster, while as given in equation (15.2) gives 

the increase of the period of oscillation on account of anharmonicity. By defining, 


F[z{t), z(0)] = \/(l + z'^{t)) - \/(l + 2^(0)) - sinh ^(l/ 2 ;^(t)) + sinh \l/2;^(0)) (16) 


Equation (15) can be finally integrated to yield (by writing t = nTo), 


F[z{t),z{0)][FD{z)/FN{z)]=-{3/2)p(^mt + M{0)t- AM{t')dt'^ (16.1) 

where z{t) thus denotes the variation of 2 when change over is made from discrete n 
to continous variable t. Using the above equation (16) we plot v{t) with t for globular 
clusters and open clusters. 

As a first step we solve equation (16) by taking as a zeroth order approximation 
Fn = 1 = Fd. In the next step we solve equation (16) again by taking Fn{z) = FN{zO{t)) 
and Fd{z) = FD{zO{t)). The resulting z{t) values are next used to calculate Fm{z) and 
Fd{z), and solve for the new z{t) by using these values of Fn{z) and Fd{z) in equation 
(16). This iterative scheme of calculation is now being extended to take care of iteration 
upto any desired level. Also the scheme can take care of the effects of mass loss if the 
rate of mass loss (dM/dt) can be specified for the star in different stages of its evolution. 


4.2 Case (b) : Tq > t(M); e' = To/T{My, e = e' 

This is a case when the star relaxes within a single oscillational time period. This is 
indeed an overdamped case. Here we take only the second term in equation (12.1) as the 
zeroth order term while the first term is a perturbation. This is easily integrated in the 
zeroth order approximation to be, 


F[zit),z{0)] = V(1 + z%t)) - v/(l + z6(0)) + sinh-i(l/03(t)) 

= —3/3 t + M{0) ^ ~ J {t')dt'^ 


sinh ^(l/ 2 ;^( 0 )) 
(17) 


Since in case (b) we have discarded the role of harmonic potential, it describes the motion 
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for a free particle under damping due to dynamical friction. Comparison of equation 
(16.1) with equation (17) shows that the factor 3/2 on the r.h.s. of equation (16.1) is 
replaced by 3 in equation (17). This is because when we deal with a particle in absence 
of the harmonic potential its variation in velocity is due to damping alone. For a particle 
in a harmonic (also anharmonic potential) the oscillation makes v to oscillate and the 
dynamical friction also tunes itself accordingly. In effect the damping is faster in case (b) 

The above equation (17) can be easily solved numerically. This is a solution which 
gives v{t) ~ u(0) exp(—t/r(m)) for v Um and [u(t)^ — u(0)^] = —3/3t for u Cm as is 
known for motion under dynamical friction, in absence of a potential. Next we define 

r{t) = f v{t')dt' (18) 

Jo 

which can be calculated from solution of v(t) as given in equation (17). We now use 
the numerical values of r(t) for the first term in equation (12.1), being thus treated as a 
perturbation. We note that the role for the first term in equation (12.1) is to provide an 
extra acceleration: 

Foit) = (19) 

We now add this extra term to equation (12.1) and integrate. This leads to 


F[z(t),ziO)] = ^/{1 + z^it))- \/(l + z6(0)) 
= -3/3 (m t + M(0) t - AM{t')dt' 


- sinh"^(l/03(t)) + sinh"\l/03(O)) 

- {l/V2amXo) Foit)dt (20) 


where the last term in equations (20) describes the effects of the gravitational field and 
is calculated from the solution of equations (17) and using the solution to numerically 
evaluate equations (18, 19) and the second term in equation (20). The r.h.s. of equation 
(20) being thus known in its entirety, the same numerical method used for solution of 
equation (17) can be repeated to find the solution of equation (20). 


5. Numerical scheme of solution 

Time development of the velocity of the star in the potential of all other stars is contained 
in the dimensionless quantity z{t) which has to be solved from equation (16) when Tq <C 
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Figure 1. (a) F[z(t),z(0)] for different z(0) values, (b) Approximations for xFch{x),xF{x), 
<&(?/) for the globular cluster case (c) Normalised gravitational field in the open cluster 
(NGC 2099) and a typical globular cluster. 

t{M) and from equations (17) and (20) for Tq t{M). It is , however, necessary to 
note that in both cases the formal schemes are the same. Essentially, we have to solve 
equations of the type. 


f{z{t),z{0j) = x{t) 


( 21 ) 


In case (a) f[z{t),z{0)] = F[z{t), z{0)][Fd{z)/Fn{z)] while in case (b) f[z(t),z(0)] = 
F[z(t),z(0)] and the x{t) are the r.h.s. of equations (16), (17) or (20) as the case may be. 
The numerical scheme of solution is thus essentially the same and is described below. 

Given any initial condition z(0) we tabulate the l.h.s. of equation (21) for different z 
and x(t) is also caculated for different t. It is then easy to find graphically the point z{t) 
at which f[z(t),z{0)] equals x(t) which is shown in figure 1. One must, however, observe 
that while evaluation of the l.h.s. of equation (21) is a simple numerical work, evaluation 
of the r.h.s. of equation (21) needs inputs on mass loss from various sources. 









14 


S. Chatterjee et al. 


6. Calculation of AM(t) 

We assume t to begin from ths star formation stage and AM {t) must contain the mass loss 
from all the sources e.g. mass burning, winds, mass ejection etc. A notable feature of the 
present theory is that the mass appears only on the r.h.s. of the equations (16), (17) and 
(20), while the dynamical variables like v appear only on the l.h.s. This greatly helps the 
numerical work and yet the evolutionary history and dynamical history of the particle are 
both accommodated in one single equation. It is clear that the mass reduction in the main 
sequence occurs due to mass burning of upto 10% of core hydrogen, converting only a 
fraction of 0.0007 of the total hydrogen mass into energy. The quantity AM{t')dt' due 
to mass burning can thus be neglected in comparison to M(Q)t, within the main sequence 
life time of the star. In this phase however, some mass is also lost due to stellar winds, 
the total mass loss due to mass burning and stellar winds is taken from Schaller (1992) 
for solar composition and the quantity AM{t')dt' can thus be computed. This goes as 
an input to the r.h.s. of equations (16, 16.1, 17, 20). Beyond the main sequence life time, 
the mass loss AM(t) becomes extremely complex for analytical calculation (Wilson 2000, 
Kudritzki and Puls 2000, Iben and Renzini 1983, Iben 1974). However the AM{t) values 
provided by Schaller (1992) incorporates all these processes for stars in various phases of 
their evolution and has been incorporated to find f* AM{t')dt' in our computation and 
thus to find the dynamical evolution of stars of different masses present in the cluster. 


7. On the question of mass segregation 

The role of dynamical friction is to slow down fast moving particles and to finally equi¬ 
librate the system so that the velocities of the different masses reach their respective 
equipartition values. This energy loss by the test mass will result in a gradual shrinking 
of the amplitude a{t) of its oscillation inside the cluster. This thus results in a segre¬ 
gation of the heavier mass towards the center of the cluster. How the mass segregation 
evolves with time can be measured by noting the variation of {r(t)) with t from the virial 
theorem. By applying the virial theorem, 2(T) = {—rdV/dr), we find, (on noting that 
2(T) = 2.i(v^) = 2.i(iu^), where averaging (• • •) is done over a cycle). 


v ^/2 


(tt^/^Gpo/t) 


l{r^) 

[a + 673(r2)3]i/2 


( 22 ) 


where (valid for ) is calculated from equations (14-16) and (r) ~ (r^)^/^. The 

variation of (r) with t for different masses M is a question of importance in studying the 
mass segregation in the cluster. 
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Figure 2. (a) Evolution of uy for different masses at different initial velocities under 
dynamical friction alone, (b) Effect of dynamical friction on the evolution of velocities 
under harmonic potential. 

8. Results and discussions 

Investigations contained in the present paper give a method to study the dynamics of a 
star in a potential due to a spherically symmetric distribution of masses, - a situation 
that arises in star clusters - uy of the star being damped by the dynamical friction with 
u_L being enhanced by the stochastic heating. We have considered the mass distribution 
in the cluster to be falling off in a Guassian fashion with the distance r from the center of 
the star cluster. Such a variation is chosen because close to the center of the cluster, the 
gravitational potential is expected to be harmonic, while the anharmonic part takes over 
only at longer distances. To calculate the dynamical friction, the Chandrasekhar type 
formula which uses King type distribution of velocities of stars has been used. However 
this modified Chandrasekhar formula being quite involved, we have approximated it by a 
simple power law type of expression (13.2), which has an excellent fit with the modified 
Chandrasekhar formula, as seen from Kolmogorov-Smirnov test. It is indeed this iden¬ 
tification of a simple formula which fits the Chandrasekhar formula very satisfactorily - 
that makes our numerical scheme very simple, fast and easy to check and yet accurate. 
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Figure 3. Evolution of velocities for a globular cluster (a), an open cluster (b) The 
captions of panel (a) applies to (b) also. 



The computational scheme is further facilitated by the fact that it is adequate to approx¬ 
imate the gravitational field TV by equation (13) and the gravitational potential (/)(r) by 
equations (13.3 - 13.6). The results that follow are based on the above inputs and implies 
that we are considering an epoch after the completion of the violent relaxation phase, 
at which an initial mass segregation described by p{r) has taken place. We have chosen 
bmin = average interparticle distance and bmax = 

The results of our calculations are found in figures 1-6, which are indeed self explana¬ 
tory. In our computations we have used the numbers given in section 4 of the paper. The 
method of solution is illustrated in in figure 1(a) where we show the plot of F[z(t), 2 :( 0 )] 
with z{t) for various initial values z(0). As a check we have tried to see the evolution 
of velocity of a test particle in absence of the gravitational potential which is essentially 
the solution of equation (17). We find that in such a case the velocity decay follows 
U||(t) = v\\{Q)exp{-t/T{m)) for uy < am and [u||(t)]^ = [u||(0)]^ - 3/3t for uy > am 
(for low t values) as is expected from the Chandrasekhar dynamical friction. This is 
illustrated in figure 2(a), where only the effect of dynamical friction is retained. The 
same scheme of solution when applied to equation (17) yield the result in figure 2(b), 
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where we find that due to oscillations in the potential well, the dynamical friction gets 
modulated in terms of instantaneous velocity and hence the decay may be slower. The 
usual mass dependences are also borne out due to the dependence of dynamical friction 
on the mass M of the projectile. This also has bearing with the fact that lower masses 
lose small amounts of their mechanical energy, while larger masses lose large amounts of 
their initial mechanical energy. This has been adequately described by equations (9-11). 
However, we cannot give at present accurate v dependences of the integrals in (10) and 
(11) and hence the exact nature of approach towards equipartition cannot be described in 
the present paper and is reserved for a future work. An approximate scheme i^contained 
in (5), which takes care of dynamical friction as well as stochastic acceleration^. It shows 
the following qualitative behaviour that the decrement in v‘^{t) is more pronounced for 
the higher mass stars. This, indeed results from a competition between the two terms 
in equation (5) which shows that the net loss of energy due to dynamical friction goes 
as m{m + M) while the stochastic acceleration - in effect a diffusion in velocity space - 
goes as m? (Chandrasekhar 1943, Spitzer and Schwarzschild 1951, Wielen 1977). For a 
more exact result in a realistic situation these effects with proper mass spectrum are to 
be considered (Chatterjee 1991, 1995). This is illustrated in figure 4. These processes 
will also describe SFQ{r,t), which we have ignored in our present approximation. It can 
indeed be incorporated as a perturbation once we understand the effects of condensation 
and evaporation from our present work. 

The change of v'^it) for different masses, enables us to make contacts with earlier 
workers (e.g. Gorti and Bhatt 1996, 1996a, who consider the dynamical effects of pro- 
tostellar clumps in gas clouds and Aarseth and Heggie 1998). Their numerical works 
essentially contain the description of the evolution process as seen in equation (5). This 
evolution of n^(t) also serves as the means to investigate the mass segregation in clusters. 
Given the value of v^, the value of (r) ^ (r^)^/^ for a star of given mass M can be 
found from equation (22), the results being represented in figures 5 and 6. In obtaining 
these results we have considered only heavier stars i.e. 1 < M/Mq < 5 i.e. m <C M. 
In obtaining figure 6 we have noted that the heavier stars evolve to their late phase in 
^ 10® years i.e. within a time scale in which very little effect is observed with respect to 
their dynamical evolution. The subsequent scenarios of their evolution being not known 
with precise models, we have considered these stars to evolve only dynamically and noted 
their condensation with time, through the parameter {r(t)). In other words figure 6 sug¬ 
gests the ’’tentative” locations within which moderately heavy stars or their remnants 
e.g. neutron stars are to be found as they evolve between 10® - 10^® years. 

We return again to figure 5 and 6, to understand the possible origin of an extended 
corona in open clusters (Nilakshi et al. 2002, Pandey et al. 2001, Durgapal and Pandey 
2001). It is seen that if we have a IMq star with ti(0) = 1.5am, it can occupy upto 11 pc 


^Why equation (5) cannot describe approach to equipartition is because the two averages used in 
calculating (dv^^/dt) and d{{Av)'^)/dt have been considered to be independent. An averaging subject to 
the constraint (8) will lead to equipartition and will be discussed in a future paper. 
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Figure 4. Change in v‘^{t) for different masses when r;(0) = 1.5(Tm- 

initially (newly born) in an open cluster, while occupy only upto 5 pc after 10® years. The 
corresponding numbers for a typical globular cluster are 21 pc at birth and 20 pc at the 
end of 10^® years. If on the other hand we consider a Q.IMq star with u(0) = 1.5am at 
birth, it will continue to live within 11 pc if it be in a typical open cluster and within 21 pc 
if it be in a typical globular cluster, because of very weak velocity loss due to feebleness of 
the dynamical friction that it experiences. It is thus possible for such a star to travel far 
away from the center of the open cluster in view the smallness of the gravitational field 
that arises in an open cluster. This may thus explain the presence of an extended corona 
in older open clusters. However, if the corona are integral part of star formation process 
itself as has been indicated by Nilakshi et al. (2002), then the dynamical evolution will 
only modify it. 

The competition of dynamical friction with stochastic heating as given in equation (5) 
is indeed qualitatively correct when dealing with M ^ m but always leads to a decrease 
in u®(M) for all masses. The correct approach should be to obtain v^{M) from equation 
(11), which slows down the evolution of u®(M) as it approaches the equipartition value. 
This approach as also equation (5) in the present paper do predict the situation that 
leads to the lighter mass (a) lying very close to the edge of the corona or (b) escaping 
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Figure 5. {r{t)) versus t for a globular cluster (a) and an open cluster (b). The captions 
of panel (a) applies to (b) also. 


the gravitational potential of the cluster. Concerning case (a) we note that this has 
been borne out in several recent works (Durgapal and Pandey 2001, Pandey et al. 2001, 
Nilakshi et al. 2002) who hnd a long tail for the corona, owing to the presence of low mass 
stars . Though the density profile p{r) determined by them can not be fitted to a Gaussian 
curve and definitive form is not yet established, the long tail in the corona, can still be 
considered to give rise to gravitational held of the type xF{x) given by equation (13.3), 
which we have already tested by considering various types of p{r)'‘s with the coefficients 
a and b of the ht being decided by the nature of the p{r). The dynamics that can follow 
from such a ht is in the process of being tested by us and in a preliminary study is seen to 
be consistent with the models concluded through observations. To estimate evaporation 
process we ought to solve the Boltzmann equation by using the appropriate scattering 
cross section for the cr(u||,u_L) an inverse square held (Chandrasekhar 1943, Spitzer and 
Schwarzschild 1951), but for a particle bound in an anharmonic potential (fir) given by 
equation (13.1). Such a method will at one stroke combine (1) dynamical friction (2) 
stochastic heating and (3) diffusion in r space in an anharmonic potential <fir). The last 
mentioned has been qualitatively addressed here but a full scale understanding (e.g. of 
the density prohle and role of evaporation in describing the density prohle) does deserve 
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Figure 6. {r{t)) versus t for a globular cluster (a) and an open cluster (b). The captions 
of panel (a) applies to (b) also. 

such a combined attack. As far as we are aware, this problem has not been tackled till 
now, to which we will return shortly along with the role of an initial mass function and 
velocity dispersion on the dynamical evolution of a cluster. 
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